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The dynamics of a typical open quantum system, namely a quantum Brownian particle in a 
harmonic potential, is studied focussing on its non-Markovian regime. Both an analytic approach 
and a stochastic wave function approach are used to describe the exact time evolution of the system. 
The border between two very different dynamical regimes, the Lindblad and non-Lindblad regimes, 
is identified and the relevant physical variables governing the passage from one regime to the other 
are singled out. The non-Markovian short time dynamics is studied in detail by looking at the mean 
energy, the squeezing, the Mandel parameter and the Wigner function of the system. 
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I. INTRODUCTION 

The description of the dynamics of quantum systems 
interacting with their surrounding is, in general, a very 
difficult task due to the complexity of the environment. 
An exact approach would require to take into account not 
only the degrees of freedom of the system of interest, but 
also those of the environment. One can, in principle, look 
for the solution of the Liouville - von Neumann equation 
for the density matrix of the total closed system. Even 
when this is possible, however, the total density matrix 
contains much more information than what we actually 
need, since one is usually interested in the time evolution 
of the reduced system only. 

A common approach to the dynamics of open quan- 
tum systems consists in deriving a master equation for 
the reduced density matrix describing the temporal be- 
havior of the open system This equation is in gen- 
eral obtained by tracing over the environmental vari- 
ables, after performing a series of approximations. Two 
of the most common ones are the Rotating Wave Ap- 
proximation (RWA) and the Born-Markov approxima- 
tion. The first one basically consists in neglecting in 
the microscopic system-reservoir interaction Hamiltonian 
the counter-rotating terms responsible for the virtual ex- 
changes of energy between system and environment. The 
second one neglects the correlations between system and 
reservoir assuming that the changes in the reservoir due 
to the interaction with the system cannot feed back into 
the system's dynamics. 

The Born-Markov approximation leads to a master 
equation which can be cast in the so called Lindblad 
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form (3, 13 • Master equations in the Lindblad form are 
characterized by the fact that the dynamical group of 
the system satisfies both the semigroup property and the 
complete positivity condition, thus ensuring the preser- 
vation of positivity of the density matrix during the time 
evolution. Moreover, it has been shown that numerical 
techniques as the Monte Carlo Wave Function method 
can always, in principle, be applied to the description of 
the dynamics, provided that the master equation is in 
the Lindblad form Q. 

In many solid state systems such as photonic band 
gap materials and quantum dots the Markov approxi- 
mation is, however, not justified p|- Similarly, the reser- 
voir interacting with a single mode cavity in atom lasers 
is strongly non-Markovian These physical systems 
therefore necessitate non-Markovian analytical or numer- 
ical approaches to their dynamics. Moreover, the non- 
Markovian features become of importance when one is 
interested in the initial temporal regime, even for Marko- 
vian systems, where the memory time of the reservoir tr 
is much smaller than the system characteristic time scale 

It is worth mentioning that during the last few years 
the interest in open quantum systems has increased 
mainly due to three reasons. On the one hand, the phe- 
nomena of decoherence and dissipation, characterizing 
the dynamics of a quantum system interacting with its 
surrounding 0, are considered nowadays the major ob- 
stacles to the realization of quantum computers and other 
quantum devices @- On the other hand recent experi- 
ments on engineering of environments Q pave the way 
to new proposals aimed at creating entanglement and 
superpositions of quantum states exploiting decoherence 
and dissipation 0, 0|. Moreover, we emphasize that 
presently the fields of quantum information and open 
quantum systems actually merge in the question to what 
extent the entangling power of the system-reservoir inter- 
action is the responsible factor for decoherence [T3. Il3| . 
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In order to understand deeply the origin, essence and 
effects of decoherence phenomena, it is of great impor- 
tance to have the tools for an exact description of the 
open system dynamics. 

In this paper we focus on an ubiquitous model of the 
theory of open quantum systems, that is the quantum 
Brownian motion model, which is also known as the 




master equation obtained with the time-convolutionless 
projection operator technique pj . We exploit an analytic 
approach based on the algebraic properties of the super- 
operators appearing in the master equation [2l|. The 
knowledge of the explicit expression of the exact analytic 
solution of this equation allows us to study in a complete 
way the dynamics and, in particular, the short time non- 
Markovian regime. The analytic results are in very good 
agreement with the numerical simulations obtained by 
means of the non-Markovian Wave Function (NMWF) 
metho d |2a. a variant of the Monte Carlo (MC) methods 
II EE . By looking at the dissipation and diffu- 

sion coefficients we single out those parameters governing 
the passage from Lindblad to non-Lindblad type regimes. 
The dynamics of the system in the Lindblad type regime 
is governed by a Lindblad type master equation, whereas 
in the non-Lindblad type regime, the master equation 
cannot be cast in the Lindblad form and the dynamics 
is dominated by virtual exchanges of energy between the 
system and the environment. 

The paper is structured as follows. In Sec. |H]we recall 
the exact time-convolutionless master equation for quan- 
tum Brownian motion and its superoperatorial solution. 
In Sec. IIIII we apply the non-Markovian Wave Function 
method to the system under scrutiny. Sections IIVI and 
[V] contain the main results of the paper. In Sec. IIVI we 
study the border between Lindblad and non-Lindblad re- 
gions, and in Sec. [V]we focus on the non-Lindblad type 
dynamics looking at the temporal evolution of the squeez- 
ing, of the Mandel parameter and of the Wigner function. 
Finally, Sec. IVII concludes the paper. 



II. 



EXACT DYNAMICS OF QUANTUM 
BROWNIAN PARTICLE 



In this section we first recall the microscopic model 
and the exact master equation for the reduced density 
matrix of a quantum Brownian particle. Then we briefly 
sketch an analytic approach to derive an exact solution 
for the master equation. 

Let us consider a harmonic oscillator of frequency 
uq surrounded by a generic bosonic environment. The 
Hamiltonian H of the total system can be written as fol- 
lows 



H = H + H E + aXE, 



(1) 



tively, and a is the dimensionless coupling constant. The 
interaction Hamiltonian considered here has a simple bi- 
linear form with position of the oscillator X and position 
environmental operator E = ^ n K n x n , where x n are the 
position operators of the environmental oscillators. For 
the sake of simplicity we have written the previous ex- 
pressions in terms of dimensionless position and momen- 
tum operators for the system oscillator. The key quantity 
governing the nature of the coupling is the spectral den- 
sity J(u>) = J2 n K nS(uJ - w„)/(2m„w n ), with m„ and uj n 
masses and frequencies of the environmental oscillators, 
respectively. 



We denote by p the density matrix of the total system. 
Under the assumptions that: (i) at t = system and en- 
vironment are uncorrelated, that is p(0) = ps(0)<8> Pe(0), 
with ps and pe density matrices of the system and the 
environment respectively; (ii) the environment is sta- 
tionary, that is [He, Ps(0)] = 0; (hi) the expectation 
value of the environmental operator E is zero, that is 
Trg {Epe(0)} = (as for example in the case of a ther- 
mal reservoir), one can derive the following master equa- 
tion E2H 



dpsjt) 
dt 



in 



A(t)(x i, y -n(t)x s P s 

r(t)(X 2 ) s + z 7 (0X s P s W). (2) 



where H = uj [P 2 + X 2 ) /2, He and aXE are the sys- 
tem, environment and interaction Hamiltonians, respec- 



We indicate with X S ' E ) and P s ( s ) the commutator (an- 
ticommutator) position and momentum superoperators 
respectively and with Hg the commutator superopera- 
tor relative to the system Hamiltonian. It is not diffi- 
cult to prove that such a master equation, obtained by 
using the time-convolutionless projection operator tech- 
nique 0, 0, |3(j , is the superoperatorial version of the 
Hu-Paz-Zhang master equation 1181 . This equation can 
be solved in generality [|l E [lHUM EJ although the 
time dependent coefficients have no obvious closed form. 
For this reason, the exact study of the non-Markovian 
behaviour of the system for strong couplings may be per- 
formed only by means of numerical methods. Also in this 
case, however, one needs, firstly, to find out the parame- 
ter regime in which a perturbation expansion to a given 
order yields reliable numerical results. This is in general 
a very difficult task. In the following we will focus on the 
dynamics of the system in the weak coupling limit, and 
in particular we will consider the case in which a trunca- 
tion of the expansion in the coupling strength a to the 
second order is physically meaningful. Under this con- 
dition the time dependent coefficients appearing in the 
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master equation can be written as follows 



A(t) 
U(t) 



k(t) cos(u)ot)cIt, 



p(r) sin(woT)dr, 



k(t) sin(wo T )^ T : 



where 



and 



r(t) — 2 / p{f) cos(loot)cIt, 
Jo 



K(r)=a 2 ({E(r),E(0)}), 



(3) 
(4) 
(5) 
(6) 

(7) 
(8) 



are the noise and dissipation kernels, respectively. For 
the case of an Ohmic reservoir spectral density with 
Lorentz-Drude cutoff 1] 



J M = — 



IT UJ 2 + UJ 2 ' 

the noise and dissipation kernels assume the form 



cu c e 



k(t) = 4a 2 kTcu 2 — 

n— — oo 



-\v-n\r 



(9) 

(10) 
(11) 



where lu c is the cutoff frequency and v n = 2irnkT denote 
the Matsubara frequencies. 

It is worth noting that, as shown in Ref. it is pos- 
sible to estimate in an easy way the order of magnitude 
of the error associated to the truncated expression of the 
coefficients. This allows to check the range of validity 
of the weak coupling approximation. The errors of the 
time dependent coefficients, up to the fourth order in the 
coupling constant, are studied in Refs. |TL fl5| - 

As we will see in the following, truncating the pertur- 
bation expansion to the second order, it is possible to find 
a closed analytic form for two time dependent coefficients 
playing a special role in the dynamics: the diffusion coef- 
ficient A(t) and the dissipation coefficient j(t). Dealing 
with a closed analytic expression of these parameters al- 
lows to gain new insight in the dynamics of the open 
system. In fact, the possibility of studying analytically 
the border between Lindblad and non-Lindblad-type dy- 
namics stems from the availability of a closed expression 
for these time dependent parameters. 

Let us now look in more detail at the form of the mas- 
ter equation J2J. First of all, we note that this master 
equation is local in time, even if non-Markovian. This 
feature is typical of all the generalized master equations 
derived by using the time-convolutionless projection op- 
erator technique P, l25j o r eq uivalent approaches such as 
the superoperatorial one |23ll3l|. 



The time dependent coefficients appearing in J5J con- 
tain all the information about the short time system- 
reservoir correlation. The coefficient r(t) gives rise to a 
time dependent renormalization of the frequency of the 
oscillator. In the weak coupling limit one can show that 
r(t) gives a negligible contribution as far as the reservoir 
cutoff frequency remains finite Q . The term proportional 
to j(t) is a classical damping term while the coefficients 
A(t) and LT(t) are diffusive terms. Averaging over the 
rapidly oscillating terms appearing in the time depen- 
dent coefficients of Eq. one gets the following secular 
approximated master equation 

dps A(f) + 7 (f) r t f f 1 

— - = [a'aps - 2ap s a* + psa'a] 

- A(t) " 7(t) [aatps - 2ctpsa + /W], (12) 



the bosonic annihilation 
(X + iP) /-v/2 and a* = 



where we have introduced 
and creation operators a — 
(X - iP) /V2. The form of Eq. JH} is similar to the 
Lindblad form, with the only difference that the coef- 
ficients appearing in the master equation are time de- 
pendent. We say that this master equation is of Lind- 
blad type when the coefficients A(t) ± j(t) are positive 
at all times j^. Note, however, that Lindblad type mas- 
ter equations, contrarily to master equations of Lindblad 
form, in general do not satisfy the semigroup property. 

In what follows we focus on the secular master equa- 
tion given by Eq. 112(1. Let us stress that the secular 
approximation invoked here does not coincide with the 
RWA commonly used to describe quantum optical sys- 
tems. Indeed, as shown in Ref. |23j], differences in what 
we may call the RWA performed before or after tracing 
over the environment do exist, and they are in principle 
measurable. The RWA performed before tracing over the 
environment consists in neglecting the counter-rotating 
terms in the microscopic Hamiltonian describing the cou- 
pling between system and environment. The RWA per- 
formed after tracing over the environment is more pre- 
cisely a secular approximation, consisting in an average 
over rapidly oscillating terms, but does not wash out the 
effect of the counter-rotating terms present in the cou- 
pling Hamiltonian (see also Ref. |24[). 

It is worth noting that there exists a class of ob- 
servables not influenced by the secular approximation 
[2ll l32| . The exact time evolution of the operators be- 
longing to such a class can be obtained by solving Eq. 
(|12fl . Examples of such observables are the mean value 
of the quantum number operator (n(t)}, hereafter called 
the heating function, and the Mandel parameter Q. 

For the Ohmic spectral density introduced in Eq. JjJ, 
the analytic expression for the dissipation coefficient "/(t), 
to second order in the coupling constant is 



7(*> 



2 2 
a u>or r 

> 2 + l 



1-e 



-U! c t 



cos(woi)— re " ct sin(ojo*) , (13) 



with r — lo c /loo. 
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As for the diffusion coefficient A(t), defined in Eq. 
a simple analytic expression is obtained only in the high 
and low temperature regimes. In Appendix A we give 
the expression for A(t) for generic T and in Appendix B 
its high T and Markovian approximations. 

The master equation Q can be exactly solved by u sing 
specific algebraic properties of the superoperators |21j . 
The solution for the density matrix of the system is de- 
rived in terms of the quantum characteristic function 
jQCF) Xt(x,p) at time t, defined through the relation 

Ps(t) = ^ / Xt(x,p) e~^ x -^dxdp. (14) 

It is worth noting that one of the advantages of this ap- 
proach is the relatively easiness in calculating the ana- 
lytic expressions for the mean values of observables of 
interest by means of the relations: 

" = H)»(|rx(»,P))_ o , 

(P"> = ■ < 15) 

In the secular approximation the QCF is found to be [2l| 

Xt (x,p) = e~ A ^ 2+ ^/\ Q (e-Wx,e- r Up) , (16) 

where xo is the QCF of the initial state of the system, 
and we defined 

x = cos(wot)x + sin(u>ot)p, 

p = — sin(wo^) 2 ' + cos(ojoOp- (17) 

The quantities Ar(t) and T(t) appearing in Eq. I|16(l 
arc defined in terms of the diffusion and dissipation co- 
efficients A(f) and *y(t) respectively (see Eqs. ifHI)) and 
(JXTJ) as follows 

r(t) = 2 / j(h)dtx, (18) 
Jo 

A r (t) = e" r W f e r ^A(tx)dtx. (19) 
Jo 

Equation I jltjjl shows that the QCF is the product of an 
exponential factor, depending on both the diffusion A(f) 
and the dissipation 7(f) coefficients, and a transformed 
initial QCF. The exponential term is responsible of en- 
ergy dissipation and it is independent of the initial state 
of the system. Information on the initial state is given 
by the second term of the product, the transformed ini- 
tial QCF. In the weak coupling limit considered here the 
asymptotic values of the diffusion and dissipation coeffi- 
cients coincide with the Markovian ones (see Appendix 
B). In this case xo [e~ r ^x,e~ r ^p\ — ► 1 for long times, 
and the system approaches, as one would expect, a ther- 
mal state at reservoir temperature, whatever the initial 



state was. In general, however, for strongly coupled sys- 
tems, the steady state could be very different from the 
thermal state. For example, in Ref. [lj (see pp. 481-483) 
it is shown that, already for a/coo = 0.25 the steady state 
solution in the low temperature regime shows squeezing 
in position. 



III. NON-MARKOVIAN WAVE FUNCTION 
SIMULATIONS 

In this section we describe how to implement the non- 
Markovian wave function (NMWF) method for the study 
of quantum Brownian motion by using the stochastic un- 
ravelling of the master equation in the doubled Hilbcrt 
space |2 HH. We use Monte Carlo (MC) methods both 
to confirm the validity of the involved analytical solution 
and to demonstrate that these methods can be used to 
study the heating dynamics of a quantum Brownian par- 
ticle in very general conditions. One might think that 
it is straightforward to apply MC methods, e.g. the 
NMWF method, once the master equation of the system 
and the corresponding jump operators are known. How- 
ever, there exist situations in which the MC simulations 
become excee ding l y heavy from the computer resource 
point of view j3J, |35j. In the following we show that, 
in our case, MC methods can be used conveniently to 
study numerically the system dynamics also in the non- 
Lindblad regime where the time dependent decay coeffi- 
cients A(i) ± 7(£) may acquire temporarily negative val- 
ues. 



A. General form of non-Markovian wave function 
method in the doubled Hilbert space 

The most general form of the master equation obtained 
from time-convolutionless projection operator technique 
reads 

|p(t) = A(t)p(t)+p(t)BHt) 

+ ^CiWpWU^t), (20) 

i 

with time-dependent linear operators A(t), B (t) , Ci (t) , 
and Di (t). The unravelling of the master equation can 
be implemented by using the method of stochastic un- 
ravelling in the doubled Hilbert space l] Ti = Hs © Hs, 
where the state of the system is described by a pair of 
stochastic state vectors 

•«>-({}§)■ <M) 

The time-evolution of 8 (t) can be described as a piece- 
wise deterministic process (PDP) Q. The deterministic 
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part of the PDP is obtained by solving the following dif- 
ferential equation 



\Ji(t)9(t)\ 
\\0(t)\\ 2 



with 



and 



F W~{ B(t) 



9(t), (22) 



(23) 



Ji(t) 



Ct (i) 

o A {t) 



(24) 



where A(t), B(t), Ci(t), and Di(t) are the operators 
appearing in Eq. lf^U|l . 

The stochastic part of the PDP is described in terms 
of jumps inducing transitions of the form 



'(*) 



'(*) 



Ci(t)0(t) 



\\Ji(t)6(t)\\ \Di(t)i>(t) 
The jump rate for channel i is given by 

U(t)o(t)\\ 2 



»(*) 



(25) 



(26) 



Finally, the solution for the reduced density matrix is 
obtained as 



P (t)= / D6D9*\cj>)(<>p\P[6,t}, 



(27) 



where P[9, t] denotes the probability density functional 
and the integration is carried out over the doubled 
Hilbert space H [IH^. 



B. Implementation of the method for QBM 

The doubled Hilbert space state vector for the quan- 
tum Brownian particle reads 

where <j) n {t) and ipn(t) are the probability amplitudes in 
the Fock state basis. 

By comparing Eq. (|20|l with the master equation (|12fl . 
the operators A(t) and B(t) in Eq. I|23l) have to be chosen 
as 

A(t) = B{t) = -iu) Q a ] a- ~ {[A(t) + 7(i)]aW 

[A(t)- 7 (t)]aa t }. (29) 
Accordingly, the operators C, and D, are 



Ci(t) = = y ggPggj qt, 

C a (i) = D 2 (t) = V|A(t)+7(*)|o (30) 



and the corresponding operators J{ , given by Eq. i|24[l be- 
come 

^(^^y^( sgn[A( V 7W]a \t 

Mt) = ^ sgn [A(t) + 7(t)] a 



o 



.(31) 



When the system dynamics and occupation of the states 
is restricted to the two lowest Fock states the equations 
resemble closely the ones used for the study of Jaynes- 
Cummings model with detuning 25]. 

The statistics of the quantum jumps is described by 
the waiting time distribution function F w (t) which rep- 
resents the probability that the next jump occurs within 
the time interval [t, t + r). F w (t), derived from the prop- 
erties of the stochastic process, reads 



F w (t) = 1 - exp 



J ° i=l,2 



t ds 



(32) 



where for channel 1 (jump up, the system absorbs a quan- 
tum of energy from the environment) 

Mt) = ^Mf;(« + 1) [l<M*)l 2 + \Mt)\ 2 ] , 
II \ ) w n=0 

(33) 

and for channel 2 (jump down, the system emits a quan- 
tum of energy into the environment) 



P2(t) = 



ii^wii 2 



J2n[\Ut)\ 2 + \Mt)\ 2 }- (34) 



When the jump occurs, the choice of the decay channel 
is made according to the factors Pi(t) and P2(t). The 
times at which the jumps occur are obtained from Eq. 
(|32|l by using the method of inversion 0. 

For very low temperatures, the non-Markovian behav- 
ior of the heating function of the quantum Brownian 
particle may occur when (n) is of the order of 1CP 10 , 
see Fig. ^ To reach such an accuracy, a MC simula- 
tion for the estimation of (n) would require more than 
10 realizations to be generated. This problem may be 
circumvented by an appropriate scaling of the time de- 
pendent coefficients A (f ) ± 7 (t) of the master equation. 
The method is based upon the following considerations. 
Let us look at the properties of the Hilbert space path 
integral solution of the stochastic process corresponding 
to the unravelling of Eq. (|20|l . The Hilbert space path 
integral representation is essentially the expansion of the 
propagator of the stochastic process T [9, t\8 , t ] in the 
number of quantum jumps 0: 

00 

T[0,t\e o ,to]= £rW[M0o,to], (35) 

JV=0 

where TV denotes the number of jumps, and are the 
iV jump contributions to the propagator. As long as in 
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FIG. 1: Dynamics of the heating function (n(t)) in the short time non-Markovian regime. For the high T, graphics (c),(f),(i), 
we have used ro = uio/KT = 0.1; for the intermediate T, graphics (b),(e),(h) we have used ro = u)o/KT = 1; for low T, graphics 
(a),(d),(g), we have used ro = ujq/KT — 100. We indicate with solid line the analytic solution and with circles the simulations 
performed using the NMWF method. In the right-lower corner of all graphics we indicate whether the asymptotic long time 
value of the heating function is null (zero T reservoir) or not. 



the time evolution period of interest there is maximally 
one jump per realization, it can be shown that, in the 
weak coupling limit and for the initial conditions used 
here, the relevant contribution to the propagator is given 
by the first two terms T'°\ T^K In this case the expec- 
tation value of an arbitrary operator O is given by 

(0)(i) - J D9D6*(<p(t)\0\i>(t)) 

{t<°> [6,t\6 ,to}+TW[6,t\6 ,t ]}. (36) 

The contribution from gives the initial expectation 
value (O) (0) plus a term which is directly proportional 
to the decay coefficients A ± 7. Since is also di- 
rectly proportional to the decay coefficients we get as a 



result that the change in the expectation value is also 
proportional to the decay coefficients 

(O)(t)-<O)(0)ocA± 7 . (37) 

Thus, to ease the numerics and still to obtain the cor- 
rect result, it is possible to speed up the decay by mul- 
tiplying the coefficients A ± 7 with some suitable factor 
/3, and to do the corresponding scaling down by divid- 
ing the calculated ensemble average by the same factor 
at the end of the simulation. For the heating function 
the validity of the scaling can be seen directly from the 
analytic solution (sec Ea. (|40() of the following section). 
The scaling allows to reduce the ensemble size for the 
estimation of the heating function from the unpractical 
10 10 to the more practical 10 4 - 10 5 . 
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IV. THE LINDBLAD - NON-LINDBLAD 
BORDER 

This section contains the main results of the paper. 
Stimulated by the recent achievement in reservoir engi- 
neering techniques, we look at the dynamics of a quantum 
Brownian particle for different classes of reservoirs. We 
single out two reservoir parameters playing a key role in 
the dynamics of the open system, i.e. its temperature T 
and the frequency cutoff lo c of its spectral density. As we 
will see in this section, by varying these two parameters 
the time evolution of the system oscillator varies from 
Lindblad type to non-Lindblad type. 

A. Heating function 

In order to illustrate the changes in the dynamics of 
the system we will focus, hrst of all, on the temporal 
behaviour of the heating function (n(t)). In the follow- 
ing section we will further investigate the non-Markovian 
time evolution by looking at the Mandel parameter, at 
the squeezing properties and at the Wigner function of 
the system. 

Having in mind Eq. I|ltj[) and using Eq. I|15|) . one gets 
the following expression for the heating function 

(n(t)) = e- r «(n(0)) + \ (V r W - l) + A r (t), (38) 

with A T (t) and T(t) defined by Eqs. (JTBJ and (fTU|) . 

The asymptotic long time behavior of the heating func- 
tion, for times t much bigger than the reservoir correla- 
tion time tr = 1/lu c , is readily obtained by using the 
Markovian stationary values for A(t) and 7(i), as given 
by Eqs. CBlJ - lfB2jl. 

(n(t)) = e- rt (n(0)) + n(uo) (l - e ~ rt ) , (39) 

with n(u>o) = [e huJ °l KT — l) . This equation gives ev- 
idence for a second characteristic time of the dynam- 
ics, namely the thermalization time tt — with 
r = a 2 Ld^r 2 1 '(r 2 + 1). The thermalization time depends 
both on the coupling strength and on the ratio r = lo c /ujq 
between the reservoir cutoff frequency and the system os- 
cillator frequency. Usually, when studying QBM, one as- 
sumes that r» 1, corresponding to a natural Markovian 
reservoir with to c — > oo. In this case the thermaliza- 
tion time is simply inversely proportional to the coupling 
strength. For an "out of resonance" engineered reser- 
voir with r <C 1, tt is notably increased and therefore 
the thermalization process is slowed down. 

As we will see in the following, there exist other two 
characteristic timescales ruling the heating process: the 
period of the system oscillator t$ = l/^o and the ther- 
mal time Tth = l/^i = l/27rfcT defined as the inverse 
of the smallest positive Matsubara frequency. In Table |U 
we summarize the definitions of the four time scales we 
have introduced up to now. In general the open system 



TABLE I: Various time scales, 

time scale name symbol explanation 

reservoir correlation tr — l/u c ui c = environment cutoff 

thermalization tt = 1/r T = a 2 toor 2 /(r 2 + 1) 

system oscillator period ts = 1/luo ujq = oscillator frequency 

thermal r t h = 1/Vi vi = Matsubara frequency 

reservoir memory T corr = tr for high T 

Tcorr = Tth for medium T 

Tcorr = TS for low T 



dynamics depends strongly on the relative value of these 
four characteristic time scales (see also |3(J). 

Let us consider now the dynamics of the heating func- 
tion for times t <C Tt ■ For simplicity we consider as initial 
condition the ground state of the system oscillator. The 
generalization to a generic initial state is however direct 
and similar conclusions hold. For times much smaller 
than the thermalization time, Eq. (|38|l can be approxi- 
mated as follows: 

(n(t)) =s / (A(t x )- 7(*i))d*i, (40) 
Jo 

where Eq. I|18l) has been used. This equation shows 
that the initial dynamics of the heating function depends 
strongly on the sign of one of the time dependent coef- 
ficients of the secular master equation (|12|l . The reason 
for the heating function to depend only on the coefficient 
A(t) — j(t) and not on A(t) + r ){t) is simply related to the 
initial condition we have assumed. Indeed, when the ini- 
tial state is the ground state of the oscillator, for times 
t <C tt, the probability of a jump up [absorbtion of a 
quantum of energy from the reservoir, see Ea. H33|l ] dom- 
inates over the probability of a jump down [emission of a 
quantum of energy into the reservoir, see Ea l|34|l ]. This 
second process, which is the signature of the quantized 
nature of the reservoir, ensures the thermalization. 

Equation (|40|l shows us that the coefficient A(t) — j(t) 
is the time derivative of the heating function. Therefore, 
if A(£) > j(t) for all times t <C tt, the heating function 
grows monotonically, whereas if there exist intervals of 
time in correspondence of which A(i) < j(t), the heating 
function decreases and eventually oscillates. Wc remind 
that, for the case considered here, whenever A(t) — j(t) > 
at all times, the master equation (|12f> is of Lindblad 
type, whilst if for some time intervals A(t) — -f(t) < it 
is a non-Lindblad type master equation. 



B. Lindblad and non Lindblad regions 

To better understand such a behavior we study in more 
details the dynamics for three different regimes of the ra- 
tio r between the reservoir cutoff frequency and the sys- 
tem oscillator frequency: r » 1, r = 1 and r«l. The 
first case corresponds to the assumption commonly done 
when dealing with a natural reservoir while the last case 
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FIG. 2: Contour plots of A(r,t) = A(r,t)/2a 2 KT for high T (a), and of the coefficient A(r,t) -7(r,t) = [A(r, t) - 7(7-, t)]/a 2 
for low T (b). In (b) we have chosen 2nr c = huj c /KT = 10. In both (a) and (b) the contour line corresponding to A(r, t) = 
(a) and A(r, t) — 7(r, t) = is indicated by a thick solid line. 



corresponds to an engineered "out of resonance" reser- 
voir. 

By using Eqs. I|13l) and (|A1|) . a straightforward but 
lengthy calculation shows that, for r 3> 1, A(t) > j(t) > 
0. In this case the master equation (|12|) is always of 
Lindblad type and the heating function is a monotoni- 
cally growing function. The three upper graphics of Fig. 
H show the time evolution of the heating function for 
r = 20 in the case of low (a), intermediate (b) and high 
(c) temperatures. 

In the case of an engineered "out of resonance" reser- 
voir, that is when r <C 1, the sign of the diffusion coeffi- 
cient is positive in the low temperature regime whilst for 
intermediate and high temperatures it assumes negative 
values for some time intervals. However, for intermedi- 
ate and low temperatures, there exist intervals of time 
in correspondence of which A(t) — -f(t) < 0. Whenever 
A(t) — j(t) is negative, the heating function decreases, 
so the overall heating process is characterized by oscil- 
lations as shown in Fig. ^ (g)-(h)-(i), where the dynam- 
ics of (n(t)) for low, intermediate and high T, respec- 
tively, and r = 0.1 is plotted for < oj c t < 10. The 
decrease in the population of the ground state of the 
system oscillator, after an initial increase due to the in- 
teraction with the reservoir, is due to the emission and 
subsequent re-absorbtion of the same quantum of energy. 
Such an event is possible since the reservoir correlation 
time tr = 1/lo c is now much longer than the period of 
oscillation t s = 1/ujq. We underline that, although the 
master equation in this case is not of Lindblad type, it 
conserves the positivity of the reduced density matrix. 
This of course does not contradict the Lindblad theorem 
since the semigroup property is clearly violated for the 
reduced system dynamics 



Finally, for r = 1, one can show numerically that 
A(t) > at all times whatever is the reservoir tempera- 
ture T. Nonetheless, for intermediate and low tempera- 
tures the time dependent coefficient A(t) — ^(t) assumes 
negative values for some intervals of times. Such a situa- 
tion leads again to an oscillatory behavior of the heating 
function as shown of Fig.QJd) and Fig. H](e). 

It is worth stressing that the non-Markovian features 
of the heating function discussed here do not depend on 
the secular approximation. Indeed, as we have mentioned 
in Sec. II, Eq. I|38|l coincides with the expression derived 
from the exact solution. Stated in another way, the ap- 
pearance of virtual exchanges of energy between system 
and reservoir, characterizing the non-Lindblad region, is 
a general feature of the non-Markovian dynamics of the 
system and it is not connected with the secular approxi- 
mation. 



C. The key parameters r and T 

The border between the Lindblad and non-Lindblad re- 
gions depends on two relevant reservoir parameters: its 
temperature T and the ratio r between its cutoff fre- 
quency and the system oscillator frequency. For high 
reservoir temperatures the quantity effectively ruling the 
dynamics is the diffusion coefficient A(t), since A(t) ^S> 
7(i). In other words, for short times and high T, dif- 
fusion is always dominant with respect to dissipation. 
Oscillatory dynamics of the heating function appears for 
r <C 1 since A(i) oscillates assuming negative values. 
Figure |21 (a) shows a contour plot of A(r, t) for high T. 
The curve defined by the equation A(r, t) — for high T, 
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(a) (b) 




FIG. 3: (a) Time evolution of the position variance for an initial squeezed state having squeezing parameter s = 0.4. The 
reservoir parameters are a — 10~ 2 , r = 0.05, and r c = u c /2wKT — 0.2 • 10~ 6 . The dotted line shows the variance of the free 
harmonic oscillator in absence of coupling to the environment. The system exhibits squeezing whenever (AX) 2 (t) < 0.5. 
(b) Time evolution of the Mandel parameter for an initial Fock state |n = 3). The reservoir parameters are the same as in (a) 



with A(r, t) given by Eq. (|B4jl . defines the Lindblad non- 
Lindblad border. From Fig. [21 (a) one can see that the 
largest value of r in correspondence of which the system 
exhibits non-Lindblad oscillatory heating is r ~ 0.27. 

For decreasing temperatures the amplitude of A(t) be- 
comes smaller and smaller. Thus, also in the presence of 
an oscillatory behavior of A(t), that is when r <C 1, for 
low temperatures the diffusion coefficient remains always 
positive. In this case, however, dissipation is not negli- 
gible with respect to diffusion anymore and their com- 
bined action is such that, for intermediate and low tem- 
peratures, the non-Lindblad dynamics appears already 
for r > 1 [see Fig|5] (b)]. Stated in another way, de- 
creasing the temperature the oscillatory behavior of the 
heating function appears for higher values of the ratio 
r, which means that the non-Lindblad region becomes 
larger. Fig. [21 shows clearly that this region, correspond- 
ing to negative values of A(t) — j(t) (A(t) for high T), is 
notably wider for low T (b) than for high T (a) . 

While in this section we have investigated the border 
between the Lindblad and non-Lindblad type regions, 
in the following section we will concentrate on the non- 
Lindblad type dynamics for two reasons. The first reason 
is that, in general, due to difficulties in dealing with non- 
Lindblad type master equations, only few studies have 
been carried out in this regime. Secondly, we have shown 
elsewhere that typical non-Lindblad dynamical features 
may be experimentally revealed in the trapped ion con- 
text with currently available technology |37| . 



V. NON-MARKOVIAN DYNAMICS OF 
NON-LINDBLAD TYPE 

In this section we focus on the dynamics of a quantum 
Brownian particle in the case of interaction with an en- 
gineered out of resonance reservoir, i.e. for r <C 1. By 
using the techniques for reservoir engineering typical of 



trapped ion systems, such region of the parameter space 
is already in the grasp of the experimentalists. Indeed, 
by slightly modifying the experimental conditions used 
in Ref. |9j, the oscillatory non-Markovian dynamics of 
the heating function can be measured [37j • 

In order to characterize completely the dynamics of 
the quantum Brownian particle in the non-Lindblad re- 
gion we look at the dynamics of the squeezing, of the 
Mandel parameter and of the Wigner function for some 
exemplary initial states. We begin with the analysis of 
the squeezing in position. By using Eq. (|15|) . it is possi- 
ble to derive the following expression for the variance of 
the dimensionless position operator X: 



(AX)? 



-r(t) 



(AX) 2 cos 2 (w t) - 
C sin(2w i)] + A r (f) 
(AX) { (t)+A F (t), 



-r(t) 



(AP)l sin 2 (w t) 



(41) 



where (AX) and (AP)„ are the initial variances of po- 
sition and momentum operators, respectively, and Cq = 
\((X Q P Q + P X ) - (X )(P )) is the initial position- 
momentum correlation function. In the last line of Eq. 
(I41JI we have defined the function (AX){(t) describing 
the time evolution of the variance of a quantum har- 
monic oscillator in absence of interaction with the envi- 
ronment. In the high temperature limit, it is not difficult 
to prove that the interaction with the environment gen- 
erally causes an increase in the variance of the position 
operator with respect to its free dynamics. Indeed, for 
times much smaller than the thermalization time t <C tt, 
Eq. (|4*T|l may be approximated as follows 



(AX) 2 ~ (AX) f (t)+/'A(t 1 
Jo 



2{AX) { {t 1 ) 1 (t 1 ), (42) 



where Eqs. I|18|l - (|19|l have been used. Having in mind 
that for high temperatures the condition A(t) 3> 7(t) 
holds, one realizes that, provided that (AX) f r ee{t) is not 
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cM = 



CM = 0.07 



cM = 0.15 




-3-2-1 12 3 

co c t = 0.24 





-3-2-1 12 3 

CM = 0.32 





-3-2-1 12 3 

CM = 0.4 




FIG. 4: Contour plots of the Wigner function at different time instants. At r = the initial state is the coherent state 
correspondent to 7?e[ao] = 1, Im[ao] = 0. The reservoir parameters are a = 10~ 2 , r — 0.05, and r c = oj c /2tiKT = 10~ 7 . 



too big (corresponding to either a very high initial squeez- 
ing in position/momentum or a very high initial position- 
momentum correlation or a mixture of these cases), the 
integral appearing in Eq. (|42|l gives always a positive 
contribution. In the short time non-Markovian regime, 
for low reservoir temperatures, situations in which the 
system-environment correlations lead to a decrease in the 
squeezing in position, compared to its free time evolution, 
may in principle occur. Such a situation would be of in- 
terest since it could be exploited to generate squeezing 
through the interaction with an artificial low tempera- 
ture reservoir. We plan to investigate further this point 
in the future. 

In Fig. (a) we plot the short time behavior of 
{AX) (t) for an initial squeezed state with squeezing 
factor s = 0.4, that is with (AX)' = s/2 = 0.2 and 
(AP)q = l/(2s) = 0.8. We remind that X and P are di- 
mensionless, therefore squeezing in position corresponds 
to (AX) < 0.5. In the figure we compare the time 
evolution of the position variance for the damped har- 
monic oscillator with the case of the isolated harmonic 
oscillator. From the figure one sees clearly the effect of 
the virtual processes which tend to decrease the position 
variance and to bring it back to its initial value. The 
overall effect of the environment tends however to wash 
out the initial squeezing. 

Let us now consider the time evolution of the Mandel 



parameter [3j| 

Q( t) = U ; ( ^ u; - 1. (43) 

This quantity gives an indication of the statistics of the 
quantized mode described by the system oscillator. For 
a Fock state Q takes its lowest value Q — — 1 while for 
a coherent state Q is equal to . Therefore, values of 
Q < indicate subPoissonian statistics, while Q = 
characterizes Poissonian statistics and Q > superPois- 
sonian statistics. Using Eq. Ijl5(l we have derived the 
time evolution of the Mandel parameter as follows 

(n(tr+e-^(n(0))[Q(0)-(n(0))} 
(n(t)) 

In Fig. |3 (b) we show the time evolution of the Mandel 
parameter for an initial Fock state \n = 3). Due to the in- 
teraction with the artificial reservoir the initial temporal 
evolution is characterized by oscillations between subPos- 
sonian and Poissonian statistics of the quantized mode. 
This behaviour may be traced back to the virtual pho- 
ton exchanges between the system and the reservoir and 
therefore is typical of the non-Markovian non-Lindblad- 
type region. Looking at Eq. 144fl . and remembering that 
T(t) > , and that (n(t)) > (n(0)), it is easy to con- 
vince oneselves that, if the initial state is Poissonian or 
superPoissonian, i.e. Q > 0, the Mandel parameter will 
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remain positive at all times. In other words, the inter- 
action with the environment never creates subPoissonian 
statistics from an initial Poissonian or superPoissonian 
statistics. This conclusion is valid for generic tempera- 
tures T, provided that the weak coupling assumption is 
satisfied. 

We now look at the the system dynamics in the non- 
Lindblad regime, considering the time evolution of the 
Wigner function of an initial coherent state |ao). Having 
in mind Eq. I|16(l and recalling that the Wigner function 
is simply the Fourier transform of the quantum charac- 
teristic function, 



W(a) 



1 



(45) 



one gets 



W t (a) 



7r[A r (i) + 1/2 



exp 



laoe-r^e-^ * - a\ 
ArW + 1/2 



(.46) 



From this equation and from Fig. one sees clearly 
that the system-reservoir interaction spreads the initial 
Wigner function. Breathing of the Wigner function, that 
is the oscillation in its spread, appears in correspondence 
of the virtual processes. This is a new dynamical fea- 
ture which is absent both in the Markovian dynamics 
of the damped harmonic oscillator and in the Lindblad- 
type non-Markovian dynamics. Indeed, in both the pre- 
vious regimes, the spread in the Wigner function simply 
increases, linearly in time in the Markovian case, and 
quadratically in time in the non-Markovian Lindblad- 
type case. We note that different breathing scenarios 
for the second moments in different regimes have been 
discussed in |39j . 

In summary, the exchanges of energy between system 
and reservoir characterizing the non-Lindblad type region 
strongly influence the dynamics of the system. In general, 
if the initial state of the oscillator possesses nonclassical 
properties (as squeezing or subPoissonian statistics), the 
interaction with the environment tends to wash out such 
properties in a time scale which is dependent, as one 
would guess, on the reservoir parameters (spectral den- 
sity and temperature). In this section we have considered 
a high T engineered 'out of resonance'(i.e. with r <C 1) 
reservoir. In this case the loss of nonclassical properties 
appears in a time scale which is smaller or equal to the 
reservoir correlation time tr = 1/lo c . Moreover, the ef- 
fect of the virtual processes, which is also important in 
this time scale, may cause oscillations between classical 
and nonclassical states, as in the case shown in Fig. [3] 
(b). It is worth noting that, in the situation here con- 
sidered, the loss of nonclassical properties, as well as the 
oscillations due to virtual processes, happen in a time 
scale which is in general much shorter than the deco- 
herence time Td cc = A^/(ePa 2 ), with d separation be- 
tween the two components of a quantum superposition 
and At = h/\j2mkT de Broglie wavelenght |4Qj]. This 



is clearly related to the high temperature condition of 
interest here, which implies ^S> 1/KT. 

In order to give a more quantitative estimate one 
should, however look at some specific physical system, 
in order to fix also the other quantities appearing in the 
definition of the two time scales. To this aim we consider 
the recent experiment with trapped ions in which the de- 
coherence of different superpositions of the vibrational 
motion of the center of mass of the ion was observed [jj . 
For this system, this experiment has shown that deco- 
herence of a superpositon of two Fock states happens on 
a time scale of the order of 30/is. Very recently we have 
proposed an experiment |3?| to observe non-Markovian 
features of the heating function in the same system and 
with the same set up used in 9]. According to our cal- 
culations, by slightly modifying the experimental param- 
eters used, one could observe oscillations in the heating 
function in a time scale of the order of 3/zs, that is one 
order of magnitude less than the decoherence time mea- 
sured in that system for a superposition of Fock states. 

In conclusion, the investigation carried out in this 
section sheds light on the short time dynamics of the 
damped harmonic oscillator, focussing in particular on 
the high T regime. Further analysis of the squeezing, of 
the Mandel parameter and of the Wigner function may 
bring new insight in the time scales governing the loss 
of nonclassical properties and in their relationship with 
the decoherence and dissipation time scales. We plan to 
explore further this aspect, with particular attention to 
the low temperature case. 



VI. CONCLUSIONS 

In this paper we have studied the short time non- 
Markovian dynamics of a quantum Brownian particle 
moving in a harmonic potential. The dynamics of this 
paradigmatic open quantum system is described by a 
non-Markovian master equation which is local in time. 
This master equation cannot be recast in the Lindblad 
form. Nevertheless, under certain conditions, the master 
equation for quantum Brownian motion is of Lindblad 
type, i.e. it has the same operatorial form of the Lind- 
blad master equation but with time dependent (instead 
of constant) positive coefficients. 

In the weak coupling limit, the relevant time depen- 
dent coefficients can be cast in a closed form. In this 
case by using the exact analytic solution in terms of the 
quantum characteristic function, we have identified the 
parameters governing the passage from Lindblad type to 
non-Lindblad type master equation. These parameters 
are the reservoir temperature T and the ratio r between 
the frequency ujq of the system oscillator and the reservoir 
cutoff frequency uj c . It is worth stressing that the weak 
coupling limit we consider in the paper is of interest also 
in the light of the engineering of reservoir experiments. In 
fact, in order to observe experimentally the key features 
of the system-reservoir interaction, e.g. the role played 
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by the entanglement between system and reservoir in the 
decoherence process, the coupling between system and 
reservoir does not have to be too strong. The stronger is 
the coupling the faster is the establishment of quantum 
correlations between the system and the environment, 
and the more difficult is the experimental observability of 
their dynamics. For this reason the techniques of reser- 
voir engineering, allowing to control both the coupling 
constant and other reservoir parameters as its spectral 
density, look very promising for investigating fundamen- 
tal issues as the quantum-classical border. 

Our analysis of the short time non-Markovian region 
shows that the Lindblad type dynamics is characterized 
by a monotone increase of the heating function, and 
therefore of the energy, of the open system. In the non- 
Lindblad type region, on the contrary, oscillations in the 
mean energy of the system clearly indicate the occur- 
rence of virtual exchanges of energy between the system 
and the reservoir. Lowering the reservoir temperature in- 
creases the probability that virtual processes take place. 

It is worth noting that whenever the master equa- 
tion for the system is of Lindblad type, it is possible 
to apply the standard MC simulation schemes and there 
exists a direct correspondence between the MC simu- 
lation method and a continuous measurement scheme 
[HHH-Hl] ■ For more general non-Markovian Monte Carlo 
methods, e.g. the NMWF we have used in this paper, an 
analogous correspondence would be of interest. There 
are indications that the Lindblad non-Lindblad border 
might be identified with the border between existence 
and non-existence of a measurement scheme interpreta- 
tion for non Markovian stochastic methods. For this rea- 
son, the study of the dependence of this border from 
parameters as the reservoir temperature and the ratio r 
might give some insight and useful hints to the research 
on this contemporary topic. 

Finally, last part of our paper deals with an analytic 
description of the short time dynamics of the quantum 
Brownian particle when virtual processes dominate. We 
have investigated in detail the temporal evolution by 
looking at the squeezing properties, at the Mandel Q pa- 
rameter and at the Wigner function. Wc have found that, 
if the system initially possesses nonclassical properties as 



squeezing in one of the quadratures or non-Poissonian 
statistics, these properties tend to be washed out due to 
the interaction with the reservoir. However oscillations 
between squeezing and non-squeezing as well as between 
sub-Poissonian and Poissonian statistics appears in con- 
nection with the virtual exchanges of energy. A further 
sign of the virtual processes is the breathing in the width 
of the Wigner function. 

Summarizing, the main result of this paper is the de- 
tailed analysis of the non-Markovian features character- 
izing the dynamics of a quantum Brownian particle, with 
special attention to the appearance of virtual processes 
for certain ranges of reservoir temperature and cutoff fre- 
quency. Due to the generality of the model here studied 
we think that our results can both contribute to fun- 
damental research on open quantum systems and, when 
applied to specific physical contexts, shed light on their 
dynamics. 
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APPENDIX A: DIFFUSION COEFFICIENT 

In order to derive a closed analytic expression for A(t) 
valid for all temperatures and for all values of the ratio r 
we integrate Eq. 110| l and use the series expansion of the 
hypergeometric function: 



A(i) 



a 2 wq 



1 



1 + r 



■| coth(7r7'o) — cot("7rr c )e " c * (rcos(wo^) — sin(iJot)) 



cos(w t) [F(-r c , t) + F(r c , t) - F{ir , t) ~ F(~ir , t)] 

TTTq 

1 r e -Vit 

— sin(woi) 

7T 



L2r (l + r L . 
+ ^r(F(-r c ,t)-F(r C7 t))}} 



((r -i)G(-r ,t) + (r + i)G(r ,i)) 



(Al) 



In this equation we have used the notations ro = 
LUo/^irkT, r c = u> c /2TrkT, 

F(x,t) ee 2 Fx (x, 1, 1 + x, e'"^) (A2) 

G(x,t) ee 2 F 1 (2, l + 2;,2 + x,e- iyit ) , (A3) 



where 2F1 (a, b, c, z) is the hypergeometric function |4l| . 



13 



APPENDIX B: MARKOVIAN AND HIGH 
TEMPERATURE LIMITS 



Inserting these expressions in Eq. IjAlfl and using the 
approximations 



In the asymptotic long time limit, the time dependent 
coefficients j(t) and A(t) tend to their stationary value 
given by 



A M = A(t 

7m = 7(* - 



)) = a 2 wq 



1 



1 + r 2 



■ coth(w /2/cT)(Bl) 
(B2) 



and the master equation l|12fl becomes the well known 
Markovian master equation for a damped quantum har- 
monic oscillator 



dps 
dt 



= - T[n(Lu ) + 1] [a*aps — 2apsa) + psa^a] 
- Tn(uj ) [aa) ps — 2a) psa + psaa)] , (B3) 



with T = a 2 uj r 2 /(l + r 2 ) and n(uj ) = (e w °/ fcT - l) . 

As far as the high temperature limit is concerned, hav- 
ing in mind the definitions of ro and r c , one sees imme- 
diately that such an approximation amounts at taking 
rn, r r 1, that is x <C 1. Under this condition one has 



F(x,t) = 2 F 1 (x,l,l + x,e- Ult ) 

~ 2 Fi (x, 1, 1, e - " 1 *) = (l-e 
G(x,t) = 2 F 1 (2,l + x,2 + x,e- yit ) 



2F 1 2,l,2,e 



-Vlt\ 



1 



cot(7rr c ) 
coth(7rro) 



_ 2kT 

) 

1 2kT 



Trr 



CJ 



one gets the Caldeira-Leggett high temperature expres- 
sion for A(i) O 



A(t) 



HT 



2a 2 kT- 



1 



{1 - e~ Uet [cos(w i) 



(1/r) sin( Wo i)]} . 



(B4) 



The other time dependent coefficient, 7(f), does not de- 
pend on temperature as one can easily see by Eq. J2J. 
We stress that, comparing Eq. (|B4|) with Eq. (|13fl . one 
notices immediately that in the high temperature regime, 
A(t)» 7 (t). 
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